J-Matrix time propagation of atomic hydrogen in attosecond fields

The J-Matrix approach for scattering is extended to the time-dependent Schrödinger equation (TDSE) for one electron atoms in external few cycle attosecond fields. To this purpose, the wave function is expanded in square integrable (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document}L2) Sturmian functions and an equation system for the transition amplitudes is established. Outside the interaction zone, boundary conditions are imposed at the border in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document}L2 function space. These boundary conditions correspond to outgoing waves (Siegert states) and minimize reflections at the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document}L2 boundary grid. Outgoing wave behaviour in the asymptotic region is achieved by employing Pollaczek functions. The method enables the treatment of light - atom interactions within arbitrary external fields. Using a partial wave decomposition, the coupled differential equation system is solved by a Runge-Kutta method. As a proof of the method ionization processes of atomic hydrogen in half and few cycle attosecond fields are examined. The electron energy spectrum is calculated and the numerical implementation will be presented. Different forms of the interaction operator are considered and the convergence behaviour is discussed. Results are compared to other studies which use independent approaches like finite difference methods. Remarkable agreement is achieved even with strong field strengths of the electromagnetic field. It is demonstrated that expanding in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document}L2 functions and imposing boundary conditions at the limit in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document}L2 function space can be an advantageous alternative to conventional propagation methods using complex absorbing potentials or complex scaling.

The availability of tailored electromagnetic fields in the attosecond range 1-3 allows examining time-dependent phenomena that don't appear in the long time limit. Carrier-envelope phase effects in ionization 4-6 and recombination 7,8 processes and the control of wave packets 9,10 with combined UV and IR pulses 11,12 can both be named as examples and even half cycle pulses in the attosecond range seem to be realizable 13 . On the theoretical side, describing ionization processes is a challanging task, because the wave function propagates without limit in space. The situation worsens with longer interaction times where the ejected particle can escape to infinity. In almost all cases, a partial wave expansion of the Schrödinger equation is used. The angular part is represented by spherical harmonics and the radial part is discretised; techniques like B-Splines 14 , finite difference methods 6,15 or the discrete variable representation 16,17 are employed. Strong fields and long interaction times lead to occupation of high partial waves. Hence, the radial position space above a given limit is disregarded or, alternatively, a limit is chosen which is so big that, throughout the whole time of propagation, it is not exceeded. Other approaches use Dirichlet to Neumann or Neumann to Dirichlet boundary conditions at some specified radial distance r = R in the spatial domain to solve the TDSE or formulate boundary conditions via a convolution integral with the Feynman Propagator of the interaction free Schrödinger equation [18][19][20][21] . A further method to deal with limited radial distances is complex scaling. There, the wave function is rotated into the complex plane, which dampens the radial part and thereby reduces the radial extension 22 . Another possibility, which can find application with finite propagation radii or with a finite dimensional function space, is the introduction of a complex potential term localized at the boundary area of the position space or the function domain. This absorbs impinging waves to prevent interfering unphysical reflections that return into the interior zone 23 . For both methods, complex scaling as well as complex absorbing potential, no knowledge of the asymptotic structure of the wave function is required. They can be formulated independent from any boundary conditions. This is advantageous in those cases, where asymptotic behaviour at large interparticle distance is unknown or only known approximately, e.g. in scattering processes involving multi-particle continua.
In time-independent scattering e.g. electron impact ionization or in perturbative treatment of laser atom interactions, L 2 functions are commonly used. A prominent example are hydrogen-like Sturmian functions, which are used in the convergent close-coupling-and J-Matrix method [24][25][26][27][28] .
A Sturmian function approach to solve the TDSE was also employed by Hamido et al. 29 and Frapiccini et al. 30 . The focus of their work is to identify time integration methods to avoid the problems associated with the stiffness of large equation systems. They present the total wave function as superpositions of Sturmian functions with time-dependent expansion coefficients and calculate multiphoton ionization of atomic H in laser fields. In these works, asymptotic boundary conditions are not imposed explicitly: The function space is chosen to be so big that numerical effects on the boundary of the L 2 function space play no role and the boundary is not overstepped.
Boundary conditions employed in the spatial domain can not be applied directly to the function domain, they differ substantially: In the L 2 approach the function space is divided into an inner part, spanned by {φ l n } N n=0 , and an outer part, spanned by {φ l n } ∞ n=N+1 . Herein l is an integer and associated with the l-th partial wave. In the outer part, short range interactions are neglected and asymptotic solutions in the space {φ l n } ∞ n=N+1 are known analytically. Boundary conditions are imposed on the complex valued expansion coefficient a l N of the function φ l N at the border of the function space. The function φ l N is delocalized in the radial domain and is not associated with a well defined radial value r = R in the radial space. In the following part we use the term boundary condition in the L 2 function sense: An algebraic condition for the expansion coefficient at the border n = N in the L 2 space to connect the inner part of the L 2 solution to the asymptotic part.
In the case of photoionization with emission of one electron or electron scattering with excitation, the asymptotic behaviour at large distances is governed by a central Coulomb force and solutions are known analytically, both in the position space and the function space. There outgoing waves e ikr with particle impulse k and position variable r are present. Then it is advantageous, to make use of the asymptotic L 2 solutions and to incorporate these directly in the boundary conditions. Because of the norm conserving character of the TDSE, in time propagation within a finite dimensional function space -and of course in position space with finite radii -there are always reflections at the numerical boundary, which result in unphysical values of the complex expansion coefficients. Then it would be desirable to make the boundary for these outgoing parts permeable and to absorb the reflected parts with the behaviour of e −ikr . Such an approach is pursued in this work: Using a time step propagation scheme, the L 2 boundary is made permeable for at least one k-value of the impinging wave at each time step. Hadley 31 applied such an approach in the spatial domain in the context of classical electrodynamics to calculate the dispersal of electromagnetic waves in beam propagation. In the spatial domain the concept of permeability is also employed by Givoli and Neta 32 , who introduce a multiplicative form of Sommerfeld's differential equation for outgoing waves, which has to be fulfilled by the solution at the boundary. This ensures permeability for a theoretically randomly large number of discrete k-values. In practical calculations the maximum number of k-values has been limited by 3 and to the one-dimensional case.
In this article, we will adopt the concept of permeability for incident waves at the numerical boundary. In contrast to the works mentioned above, we formulate appropriate outgoing wave boundary conditions not in the position space but at the border of the L 2 function space. These time-dependent L 2 boundary conditions are based on the analytically known asymptotic solutions of the J-Matrix scattering theory and ensure permeability for a specified k-value for each time step and each partial wave. This way, at each time step the boundary can be made permeable for several different k-values. The flexibility of the method is further enhanced by adjusting the variable k in every time step. As such, the method is similar to the time-independent J-Matrix approach, whereby an expansion of the wave function in Sturmian functions is adjusted at the border of the L 2 function space to fit the expansion coefficients of the inner domain to those of the outer domain. The objectives of this article is three fold: • Establishing an equation system for the one electron TDSE based on the mathematical foundations of the J-Matrix scattering theory and incorporating time-dependent asymptotic boundary conditions. • Application to the photoionization of H in attosecond fields. This includes the numerical implementation and as a proof of the method comparing the results with values of other independent approaches. For this the H atom is exposed to few and half cycle laser fields. • To assess how well Sturmian functions can describe the ionization dynamics in exterior fields whose strength is comparable with that of interior atomic fields.
To this purpose, we will first recapitulate the connection between Sturmian functions and Coulomb functions and then set up an equation system for the expansion coefficients. Subsequently, we derive the boundary conditions in the L 2 function space. Then we explicate the numerical realization and present the results. We first examine the interaction of H with a weak few cycle attosecond pulse and then the strong field half cycle case. Next the ionization dynamics in half cycle pulses is discussed and then the interaction of H with two delayed attosecond pulses is illustrated. Finally we mention possible applications for future work. In this article the atomic units used are ( e = ℏ = m e = 1).

Theoretical methods
In time-independent collision theory, it was demonstrated in a series of papers [33][34][35][36] how Sturmian functions can be efficiently used not only to describe the bound state part but also the continuum part of the collision system. It is therefore opportune, to transfer the concept of Sturmian functions so successfully used in time-independent scattering methods to time-dependent processes involving ionization. For example in the J-Matrix theory the L 2 function space is divided in an interior and exterior asymptotic part. Coulomb-Sturmians are used to present www.nature.com/scientificreports/ both parts on an equal footing, while at the boundary between the two parts the wave function is adjusted by the L 2 analogon of the R-Matrix 28 .
To begin with, Sturmian functions -which are the basis of the J-Matrix method -will be presented in short and then the algebraic description of the Schrödinger equation will be developed. This results in a system of differential equations of first order in time, which will serve as the basis for specifying boundary conditions and the numerical implementation.

System of equations for the time propagation.
The three-dimensional TDSE for an atom with core charge z is: Herein V ( r, t) is an exterior electric field in dipole approximation, either given in length, velocity or acceleration representation.
Herein E symbolizes the electric field strength, A the vector potential and the polarization ǫ of the electric pulse is given through � ǫ = � E/ E , α is the excursion amplitude of the electron experienced by the field E. For (1) a partial wave decomposition is carried out and the wave function is expanded in the radial part in Coulomb-Sturmian functions φ l n (ξ r) and spherical harmonics Y lm in the angular part: The expansion is cutoff at an upper value N which characterizes the L 2 boundary. The time dependence is contained in the expansion coefficients a lm n (t) ; the Sturmian functions are defined as follows: [34][35][36] These are regular at the origin and consist of Laguerre polynomials L 2l+1 n , powers in r and decaying exponentials. ξ -within certain limits -is a freely selectable scaling parameter in radial direction. The function's maximum is found to be roughly at 2n ξ . As a result, by varying ξ , various areas in position space can be sampled. The functions φ l n form a complete set, are not orthogonal, but tridiagonal, whereby T l n n ′ = dr φ l n φ l n ′ only differs from 0 when n ′ = n−1, n, n+1 . The integrals can be analytically evaluated: T l n n+1 = T l n+1 n = −(n + 1) 2l+2 /ξ and T l n n = 2(n + l + 1)(n + 1) 2l+1 /ξ , with the Pochhammer symbol (a) b = Ŵ(a + b)/ Ŵ(a) . The tridiagonality of the matrix elements also applies to the radial Hamilton operator with angular momentum l: 36 The expansion coefficients a lm n will be determined using the Galerkin method on the one electron basis Y lm φ l n , whereby n lies in the interval [0, N]. Projection on this basis leads to a coupled system of differential equations of first order in a lm n , thereby minimizing the deviation in the weighted mean for the expansion-coefficients: Special attention has to be paid for the boundary N, which should be high enough to enclose the major part of the interaction term V ( r, t) . For n ≥ N it is assumed, that the net photon absorption tends to ∼ 0 , meaning that the electric field can be neglected. In addition there the time derivative of a lm N+1 is disregarded, leading to the following approximation at the boundary n = N: www.nature.com/scientificreports/ Solving this equation according to i ∂ ∂ t a lm N gives: Insertion into (6) with index N − 1 leads to: The matrix elements T l n n ′ are time-independent, equation (6) with radial indices n ∈ [0, N − 2] together with (9) with n = N − 1 defines N equations for the N + 1 unknown coefficients a lm n for each partial wave. The Galerkin method is independent of boundary conditions, up to now, no boundary conditions had entered into equations (6) and (9). A single degree of freedom remains, which can be shaped so as to create outgoing wave conditions at the boundary limit n = N . This condition and the calculation of a lm N will be derived below.
Boundary conditions. If the potential term V in (1) is ignored, the solutions are given by Coulomb functions with energy E = k 2 /2 . These exhibit the asymptotic behaviour e +i(kr+ 1 k ln 2kr+ l π 2 ) and e −i(kr+ 1 and form a fundamental system of H 0 . The corresponding Coulomb solutions + l and − l can be expanded in the infinite dimensional Sturmian basis φ l n and form the foundation of the J-Matrix method: 28,[34][35][36] The expansion coefficients Q +l n consist of the Coulomb spectral-function +l 0 and the Pollaczek functions q +l n , defined by and γ , χ and x are defined by The second solution of the fundamental system (� + l , � − l ) arises via substitution k → −k . It needs to be mentioned that functions + l and − l only merge with corresponding Coulomb functions in the asymptotic limit r → ∞ when n-values are large, because with the L 2 basis φ l n the irregular behaviour ∼ 1/r l+1 at the origin cannot be reconstructed.
The various forms of the dipole operator show different behaviours for r → ∞ . The acceleration operator scales like 1/r 2 for r → ∞ , whereas the length and velocity form grow unhindered. In order to nevertheless formulate boundary conditions in the asymptotic function space for these cases, it is hypothesized that for large n-values the effect of the dipole operator can be neglected -this is equivalent that absorption and emission of photons occurs only in a restricted area around the nucleus. This approximation is strictly true only for the acceleration form, where only a restricted domain around r = 0 contributes in the interaction, but is also adopted here for length and velocity representation. With this approximation and assuming an unlimited propagation space, the coefficients a lm n in (3) have to show outgoing wave behaviour for large n and each partial wave and are thus proportional to the Pollaczek-functions Q +l n appearing in (10a): For sufficient large indices n the time evolution is entirely governed by the free field operator H 0 and thus the quotient of the expansion coefficients a lm n−1 /a lm n in consecutive n − 1, n can be regarded as approximately equal. Then the following two equations at the boundary can be set up for the quotient of two consecutive coefficients: Herein, for pure outgoing behaviour, k is real and > 0 . Due to reflections at the numerical boundary, k can take on complex values: The case Re(k) < 0 corresponds to the unphysical event where the impinging wave is scattered back. So the procedure to avoid reflections and guaranteeing outgoing waves consists of setting Re(k) = 0 when Re(k) < 0 in each time step and each partial wave. With the k-value thus determined and adjusted if Re(k) < 0 , the expansion coefficient a lm N can then be computed via the equation (12b): In the high energy limit k → ∞ the Pollaczek-function can be analytically evaluated. This yields the limit condition: Equations (13) and (14) are solved at each time step and for each partial wave. Thus k can vary and can be adjusted in each time step for each l. Using this procedure the probability of the electron to be in the interior zone n < N of the L 2 function space decreases. This is also evident in the numerical simulations, where the test cases show a decreasing probability, dependent on the basis size and the scaling parameter involved.

Role of gauge.
All three forms of the transitional operator (2) are equivalent given an infinite dimensional function space ( N→∞ ) and are connected via a gauge transformation. Because the operators are vectors, they only connect angular momenta differing by ±1 and states of different parity. Differences occur in their action in configuration and function space: The acceleration form corresponds to the Kramer-Henneberger's reference system and results when a Taylor expansion of 1 |� r+� α| for a small quiver amplitude α is taken. Thereby regions close to the core are sampled and for r → ∞ the asymptotic behaviour is characterized by 1 r 2 . In contrast, the length and velocity forms grow without limit for r → ∞.
Like the matrix elements for H 0 , the elements for the length and velocity form of the dipole operator show a narrow banded structure, while those for the acceleration form are not banded. The operators V L and V V connect only adjacent indices n, while V A leads to transitions to all Sturmians with n ′ ≥ n : With the length form, only orbitals with n ′ ∈ [−3, +3] can be reached. Per time step the wave function spreads outward by three φ l n units. With the velocity form, the wave function propagates in time by two functions in the L 2 function space.
Despite the fact that in position space the acceleration form is well localized at the origin, it spreads much more rapidly in the Sturmian function space. The behaviour is determined via the matrix element ∞ 0 drφ l n 1 r 2 φ l+1 n ′ . This is = 0 for all n ′ ≥ n , i.e. in time propagation the whole φ l n space will already be occupied in the first time step, although it will decrease as n ′ -values rise. The acceleration form is nonlocal in the indices n, n ′ and therefore less suited for time propagation. Evidence for this is also supported by numerical tests: In order to achieve www.nature.com/scientificreports/ convergence for the acceleration form, N has to be chosen much larger than in the length or velocity case with enormous costs on computer ressources. Summing up: The length form and the velocity form are suitable for time propagation and provide comparable results to circa three digits. The acceleration form, however, is less suitable for the numerical time propagation within a Sturmian function approach.
Observables. After turning off the external field, relaxation occurs and the system can be described as a superposition of Coulomb waves. All information about ionization processes is contained in the expansion coefficients a lm n (t) . Physical observables like energy spectra are obtained by projecting �(r, T f ) on Coulomb functions − k . These are normalized to incoming wave boundary conditions and correspond to an ionization experiment in which the energy of the photoelectron is measured after the ionization process.
The Coulomb wave − � k can be expanded in spherical harmonics in the impulse and position directions and an energy-dependent radial function L E .
Herein σ l (k) is the Coulomb scattering phase, given as arg [Ŵ(l + 1 − i z k )] , and l E corresponds up to a phase factor to the regular Coulomb function, normalized to an energy delta function. 37 l E can, as is shown in the J-Matrix theory, be expanded according to the L 2 functions φ l n : 36 with +l 0 and x given by (10) and p l n is a Pollaczek-polynomial 36 . The function b k is a measure how much a Coulomb wave with energy E and impuls k contributes to the solution �( r, T f ) . The impulse distribution of the photoionized electrons can be further separated into an energy part c lm E with E = k 2 /2 and an angle part. The angle part of the electron impulse is characterized by spherical harmonics and expanded as following: Insertion in (17) results in the following representation: Integration over the angle part d k and projection with �Y lm l E | gives: The angle integrated cross section dP/dE, differential in energy, is defined as: The coefficients c lm E can be determined analytically using the tridiagonal structure of the matrix elements of φ l n . The result is

Numerical realization
The equation system (6) and (9), together with the boundary condition (13) and (14), represents a differential equation system of first order for the expansion coefficients a lm n , which is to be solved numerically by a Runge-Kutta procedure. Written in vectorial notation: www.nature.com/scientificreports/ The matrix T is tridiagonal, time-independent and consists of the coefficients from the right side of (6) and (7). A transformation y(t) = Ta(t) leads to the standard form: The Hamilton matrix H in the basis of φ l n is block diagonal in l, the Hamilton matrix without laser interaction is tridiagonal and the dipole operator is pentadiagonal (depending on whether the length or velocity form is used). The time propagation is implemented with a Runge-Kutta method of 8 th order with constant step width, depending on the dimension of the equation system. At each time step a matrix vector multiplication T −1 y has to be carried out, and the resulting vector is then multiplied with the Hamilton Matrix. Due to the sparse occupation of H, only a few calculation steps are involved. In order to achieve outgoing waves at the limit of the L 2 function space it is necessary, at every time step, to determine the value a lm N pursuant to (14). To this purpose, the complex value of k is determined in accordance with (13). The Pollaczek functions q +l n are complex-valued and are calculated by means of the Gauss continued fraction. Due to the partial wave expansion, the corresponding k-values satisfying (12a) can be adjusted specifically for each l. So at each time step, the boundary can be made permeable for several k-values, thus enhancing the flexibilty of the boundary approach (11). The k-values are determined by a Newton procedure to determine a zero of (13). Because the k-values only change slowly from one time step to the next, the k-values of the previous step are used as a starting point for the Newton procedure, resulting in fast convergence after 3 -4 steps. As can be seen in (12a) for determining the k-values, there is no coupling between different l. This makes it possible in the computer code to do the calculation in parallel on a multiprocessor system using a multithreading library.
The differential system (6) is stiff, resulting in decreasing time steps with increasing basis size. This is because the positive energy spectrum is not bounded. Increasing the basis size results in approximating more and more higher energy states, which oscillate strongly dependent on the energy. Thus for longer interaction times and huge basis sizes , a Runge-Kutta method is not the favourite method, in this case Lanczos procedures or a method proposed by Fatunla 30,38 can be considered as alternatives. A further approch to overcome stiffness and to extend the method into the femtosecond range would be to transform to a spectral representation by diagonalizing a large Hamilton matrix and to disregard the high energy eigenvectors. In the simulations presented in the following part the time steps are adjusted properly and differ in the specific calculations dependent on the maximum number of functions used in the expansion (3).
All the calculations were carried out either in with the velocity or the length form of the dipole operator. The results using length or velocity form differ only minimally from one another and agree within the first 2 -3 digits. This is shown in Table 1 where the angle-integrated ionization probabilities at different energy values of the continuum electron are presented; length form shown in 2 nd column, velocity form in 3 rd column. The calculation is based on an interaction with E(t) = E 0 sin(ωt) with E 0 = 1, ω = 2 and a duration T = 2π ω = π ; it was calculated with 16 partial waves, each with 75 radial functions. For the results presented below a software package (Apfloat) was used for the numerical implementation. It permits calculations of arbitrary precision. Typically, 16-25 digits were used. In order to achieve reliable, i.e. convergent results, the maximal number of partial waves, the number of radial functions N and also the scaling parameter ξ were varied. The results shown in the following section can be considered to be convergent with respect to these parameters.

Results and discussion
Apart from the stability of the computed values with regard to variation of the parameters like basis size, numerically predicted physical observables also have to be consistent with those of other work. This is pursued in this section where the quality of the L 2 expansion is demonstrated. First, two cases are examined: 1. In the weak field perturbative regime, the method is applied to a single attosecond pulse and results are compared with those of Della Picca 39 et al. who use a differencing scheme on a spatial grid.  www.nature.com/scientificreports/ 2. In the strong field case, the ionization of H in a half cycle pulse is investigated and compared with results of Duchateau et. al. 40 who apply a B-Spline approach. The study is continued examining the dynamics of H in half cycle pulses and in a 12 cycle pulse. The results are compared with the First Magnus Approximation of Dimitrovski et al. 41 .
As a further application of the method, interference phenomena in two delayed attosecond pulses are analyzed.

Comparison with B-spline and finite difference methods.
To test the quality of the expansion in the perturbative regime, the same laser parameters used by Della Picca et. al. 39 are chosen: Herein ω = 1.71 , the number of cycles = 7, the total duration T f = 25.72 , field strength E 0 = 0.05 with linear polarization in z-direction. Fig. 1 shows the ionization spectrum of H with the ground state as the initial state. For the calculation, the velocity form was used and the convergence with respect to variation of the scaling parameter ξ and the basis size was tested. The electric field amplitude is considerably smaller than the inner atomic field strength. Thus, the one photon absorption dominates in the ionization spectrum, which can be clearly seen around E = 1.21 , where a maximum appears. At 2.9 a.u., the second peak turns out considerably smaller and corresponds to a two photon absorption. Due to the duration of 25.72 a.u., the laser pulse is of course not monochromatic, but spectrally distributed across the interval ω = 4π T f ≈ 0.5 . This also explains why the peaks at 1.21 and 2.91 are relatively wide. The significant structures in the area of the maximum at E= 1.2 are reproduced very well even with a basis size of 4 partial waves and 60 radial functions. In contrast, the structures in the area above 3.3 are smaller than at the maximum at E = 1.2 by a factor ∼ 10 −7 . In order to resolve this area, it was necessary to adopt equations with eight partial waves and 240 radial functions in the L 2 expansion. To give an impression about the convergence behaviour, calculations with different basis sizes for the cutoff parameter N in the expansion (3) are included. In all calculations 8 partial waves were employed. For N = 60 -thin green curve -the approximation is poor and only at the maximum at E ∼ 1.3 the result is reliable. For higher N (N = 120, 180, 240) the convergence behaviour improves, especially in the higher energy part above 1.8.
For comparison, results from [39] were integrated and shown by the red dotted solid curve . The minor differences above 3.5 a.u. can eventually be explained by the limited basis size and correspond to the high energy part of the wave packet which propagates outside in radial direction very quickly and thus reaches larger radial distances. The agreement -aside from the very weakly pronounced structure above 3.5 a.u -is remarkable and shows that a hydrogen like Sturmian basis can very well approximate the ionization dynamics -at least in the perturbative regime. It is, of course, of particular interest just how accurate the processes can be described when the strength of the exterior field is of the same magnitude as the inner atomic field strength. To this purpose, calculations were carried out with the H atom exposed for a half cycle to an intense field with amplitude E 0 = 1 . To ensure comparability, the same field parameters as in Duchateau et.al. 40 were taken: To sum up: From the time-independent J-Matrix method it is to be expected that an expansion in L 2 Sturmians, normally used for bound states, should be well suited to describe ionization processes in weak external fields -this is confirmed in Fig. 1. Furthermore, Fig. 2 gives convincing results that the L 2 expansion can also be applied to strong field processes in which the division of the system into 'atom' and 'outer field' no longer seems appropriate.
Building on these results, the following section examines the ionization dynamics of H in 0.5, 1, 1.5 and 12 cycle fields.

Dynamics of H ionization in half and few cycle pulses.
In this part, the effects of a half cycle pulse (HCP) on photoelectron spectra are considered. HCPs are a field of actual research, theoretically as well as experimentally 13,[42][43][44][45] . A HCP is an idealization and does not exist in reality; the time integral over the electromagnetic signal has to be zero. A HCP can be approximated by an unipolar one cycle pulse where the electric field in the first half cycle ist strongly peaked while the second half cycle is strongly prolonged with a much lower amplitude.
To begin with, the H atom, initially in its ground state, is exposed to one HCP with linear polarization of the form: Fig. 3 shows the ionization spectrum after passage of a HCP with duration T f = 0.5 2π ω for different values of the electric field amplitude: E 0 = 0.1, 0.5, 1.0, 1.5, 2.0 . These values comprise intensities which can be treated via perturbation theory, and values, which are comparable to, or larger than the inner atomic field strength. From Fig. 3 it becomes evident that for larger values of the field amplitude E 0 a significantly pronounced maximum appears that begins for E 0 > 1.
This behaviour stands in excellent congruence with the short time approximation, i.e. First Magnus Approximation, developed by Dimitrovski et al. 41 . In this approach, the transition probability for short interaction times (compared to the atomic orbital time of electrons) is expressed through: The transition probability is independent of the exact form of the electric pulse and only dependent on the time integral over the elctric field � q = t 0 dt ′ � E(t ′ ) , which describes a momentum transfer averaged over time. The same matrix element also appears in atomic scattering processes in the first Born approximation. In this case, q designates the momentum transfer of the incident particle on the target atom. The final state wave function can (28)  www.nature.com/scientificreports/ be approximated for larger kinetic energies and, therefore, for large momentum transfers q, by plain waves. For the matrix element this results in 41 : It is this very behaviour that Fig. 3 illustrates, although here it is calculated by solving the three-dimensional Schrödinger equation with an L 2 approach. According to (30), the maximum ionization probability is E = q 2 2 − |E i | , where E i denotes the ground state energy of the H atom. For q < 1 this maximum is in the negative energy range, which doesn't contribute to the ionization spectrum. For q > 1 the maximum shifts further outward into the positive, and therefore measurable, energy range. The correspondence of the 'exact' solution of the TDSE, and the results as predicted by the First Magnus Approximation, is excellent. In Fig. 3, the black curve with field amplitude E 0 =1 corresponds to the q = 1 case of the short time Magnus approximation, these are shown as red dots and are nearly congruent. Figure 4 shows the electron density along the z-axis after passage of one HCP and two HCPs, each with linear polarization pointing along the z-direction and again with the H atom initially in its ground state. In the one HCP case the density is strongly located at negative z-values with a maximum at z ∼ 3.0 . When one observes a pulse consisting of a full oscillation cycle, the electric field in the second half cycle reverses its sign: the electron is slowed down and electrons with low energies are accelerated back toward the nucleus; there they can eventually recombine to bound states. This behaviour is clearly evident in the full cycle case shown by the dashed curve in Fig. 4. The electron density is shifted to smaller negative z-values. To sum up, the electron density shows, that,  www.nature.com/scientificreports/ even after a symmetrical pulse (here, after a whole cycle) and a spatial isotropic initial state (here 1s ground state), the system is asymmetrically localized along the electric field direction. In Fig. 5, the energy spectrum after interaction with one, two and three HCPs is presented, each for different field amplitudes E 0 . In the E 0 = 0.1 case ( upper panel) and after passage of one HCP, electrons are mainly emitted with energies ∼ 0 and the probabilty drops down during the second cycle, while for electron energies > 0.3 the probabilty increases slightly -see the black curve. The case E 0 = 0.1 corresponds to the perturbative regime, the electric field in the second half cycle is so weak, that it can only accelerate back the low energy electrons freed in the first half cycle. The situation changes with increasing field strength, as can be seen from the mid and bottom panel of Fig. 5. Here the field strength is 'high enough' in the second half cycle to effect also the higher energy www.nature.com/scientificreports/ electrons, emitted in the first half cycle. In the very strong field case with E 0 = 1.5 ( bottom panel) during the second half cycle the electron density increases around threshold and decreases for higher energies. High energy electrons freed in the first half cycle are decelerated and return back to the nucleus, but a considerable fraction has still positive energy, even after deceleration. If the number of cycles is increased, the ionization spectrum shows the characteristic peak structure of abovethreshold ionization corresponding to the long time limit: In Fig. 6, the ionization spectrum for an interaction time of 12 whole oscillation cycles with ω = 2 and E 0 = 1 is presented, this time on a logarithmic scale. The spectral width is reduced in comparison to one and two HCPs; the pulse is significantly more monochromatic. This effects the spectral energy distribution: Above the threshold peaks appear; their magnitudes decline with increasing order. They correspond to absorption of photons of frequency ω = 2 : The first at ∼ 1.5 and the other at energies E = nω − E i , where E i is the ground state energy and n specifies the number of photons absorbed.
Summarizing the results from Fig. 3, 4, 5 and 6: • The ionization spectrum in the one half cycle case shows excellent agreement with the First Magnus Approximation calculation of Dimitrovski et. al. 41 . • In the case of one whole interaction cycle, recombination occurs in the second half wave, which leads to reduced ionization, but the ionization outweighs the recombination. • Increasing the number of field cycles, the charcteristic above-threshold ionization structure -with peaks for integral multiples of the oscillation frequency -are successfully reproduced.
The methods developed in this work can be extended to more complex tailored fields, for example delayed attosecond pulses. An example of this will be examined below.
Delayed attosecond trains. In the case of several time-delayed attosecond pulses, interference patterns in the ionization spectrum occur, which scale with the time delay between the individual pulses. In this part, the ionization spectrum for interaction of two delayed pulses of the form is investigated. Here T D describes the time delay, ) and N 1 , N 2 characterize the number of cycles for pulse 1 and 2. The individual pulse durations are sufficiently short, so that electrons are emitted in a wider energy range. Assuming that the electrons are released independently in each pulse -implying that the first wave packet is not affected by the second pulse-and that they overlap in energy, the interference pattern can be understood by modelling the electron wave packets by plane waves. During the first pulse, a wave packet with energy E is created at time t ′ which develops to time t according to a 1 (E)e −iE(t−t ′ ) with amplitude a 1 . At creation time t ′ the system -initially in its ground state with energy E i -contributes a phase factor e −iE i t ′ to the wave packet. Adding this phase results in: The time delay is sufficiently short, so that the phase accumulation due to the Coulomb force in [T 1 ,T 1 + T D ] can be disregarded. In the second pulse, the electron is released at t 2 = T 1 + T D + t ′ , and the phase of the ground state develops according e −iE i (T 1 +T D +t ′ ) , giving : www.nature.com/scientificreports/ The two wave packets overlap in energy. Superposing A 1 and A 2 and evaluating the squared modulus gives for the ionization spectrum: The maxima and minima are determined through the argument of the cosine function and are independent of the electric field strength. Figure 7 shows the ionization spectrum after interaction with two delayed attosecond pulses for two different field amplitudes. The electric field has the form of (31) with ω = 2 and E 0 = 1 (black curve) and E 0 = 2 (red curve). Each pulse consists of two cycles and the time delay is a whole oscillation period, T D = π = 2π ω . By way of illustration, pursuant to (34), the maxima are incorporated and represented by dashed vertical lines. The gap is given by 2π/(T 1 + T D ) and provides an excellent correspondence of the modulation pattern with the predictions from equation (34). When E 0 = 2 , the ionization probability is greatly enhanced. For both field strengths, maxima and minima are shown at the same position and fit surprisingly well with (34). The agreement of the interference structure with this plane wave approximation suggests the follwoing picture -at least in the strong field limit: • The Coulomb interaction plays only a minor role; the interference structure can be described by plane waves.
• The wave packet produced by the first pulse is hardly affected by the second pulse (single photon regime). The oscillations in the energy spectrum can be seen as amplitudes and phase dependent overlays of the individual wave packets. • Vice versa, the time delay of the pulses can be reconstructed from the interference patterns, while the height of the peaks scale with the field strength of the pulse.

Conclusion
The aim of this article was to realize a TDSE solution with asymptotic boundary conditions based on L 2 functions which minimize reflections. These boundary conditions were formulated as outgoing waves based on analytical functions of the J-Matrix scattering theory. The method presented herein is competitive with other techniques like B-Spline or finite difference procedures and gives consistent results.
• The L 2 approach makes explicit use of the known analytic asymptotic solutions, which are fitted to the coupled differential equation system of the TDSE. The numerical implementation is significantly more efficient: at every time step the Hamilton matrix is only sparsely occupied. The affected matrix elements can be calculated a priori and standard procedures can be implemented for the time propagation. • The determination of an optimal value for minimizing reflections can be realized adaptively in every time step, but does require a higher effort to calculate the Pollaczek functions. • Physically relevant observables can be extracted by projecting the time dependent wave function onto analytically known Coulomb functions of the J-Matrix scattering theory. (32b) (34) E min = E i + π T 1 + T D (2n + 1), E max = E i + π T 1 + T D (2n), n = 0, 1, 2 . . .

Figure 7.
Ionization spektrum of H after irradiation with two delayed pulses according to (31). The electric field of the pulse is depicted in the inlay with N 1 = N 2 = 2, ω = 2, T D = 2π/ω . The curves correspond to different field amplitudes: black curve: E 0 = 1 ; red curve: E 0 = 2 . The vertical lines illustrate the positions of the maxima predicted by (34). www.nature.com/scientificreports/ A priori, it is to be expected that, with an expansion according to L 2 functions, this basis will describe the underlying phenomena well such as multiphoton ionization or above-treshold ionization in the weak field area which is also accessable via perturbation theory. A further aim of this article is to explicate the quality of L 2 expansion in areas in which a perturbation theoretical treatment no longer applies, i.e. in areas, where the strength of external fields become similar to or stronger than the inner atomic fields. The results can be summed up as follows: • The method presented in this article allows one to describe ionization phenomena which are no longer accessible according to perturbation theory. The findings compare well with other independent methods. • For short times, there is an excellent congruence with other independent approaches such as the B-spline and short time Magnus approximation. • Interference structures in the ionization spectrum, based on delayed attosecond pulses, are successfully reproduced.
Only a small proportion of processes in attosecond fields were examined in this article. The method discussed here is well suited to describing other interesting phenomena: • Description of the energy spectra, when phase effects of the external few cycle field are considered (Carrier Envelope Phase). • Description of atomic systems within the single active electron approach. There the action with the core is incorporated through an effective potential, which in many cases contains a screened Coulomb potential and sums of exponentials combined with powers in the radial coordinate. Such potentials can be represented by Sturmian functions and can be implemented in the method presented herein. • The method presented here can principially used for more exotic tailored external fields, e.g. vortex laser fields. • The simulation of scattering effects in a diatomic molecule ion in the Born Oppenheimer approximation. This can be considered by making allowance for a further potential term in the TDSE.

Data availability
The numerical datasets generated during the current study are available from the corresponding author upon reasonable request.